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Recent advances in Random Matrix Theory enable one to determine the pseudoscalar decay 
constant from the response of eigenmodes of quenched fermions to an imaginary isospin chemical 
potential. We perform a pilot test of this idea, from simulations with two flavors of dynamical 
overlap fermions. 



I. INTRODUCTION 



o 
o 

The leading-order chiral effective Lagrangian 

^ i C cS = ^Tr(^CWt) - ~Tt[M(U + C7+)] (1) 

is characterized by two coupling constants, F and S. While they are fundamental from the point of the chiral 
Lagrangian, their values may be computed from the parent theory, QCD. This is a long-standing goal of lattice 
simulations. 

F and £ are not measured directly in experiment: the values of the decay constant and condensate for nonzero 
quark mass are computed from the chiral Lagrangian through familiar functions of F and S, plus contributions from 
i-C . higher order terms in the chiral Lagrangian. Lattice simulations regularly present results for the mass-dependent 
quantities. Since it is not so easy to do simulations directly at zero quark mass (or at the physical up or down quark 
masses) , the analysis of simulations typically involves fits of lattice data at several values of the quark mass to chiral 
^ . perturbation theory formulas. F and X are parameters determined as part of the fit. 

| Many techniques for extracting the mass dependent parameters, or F and £ themselves, have been devised. Most 
t ■ of them require studying the asymptotic behavior of correlation functions, and hence, need larger simulation volumes. 
However, in the epsilon regime (m^L <C 1 but AL 3> 1 for any other nonperturbative energy scale A), properties 
of the chiral Lagrangian are encoded in the spectrum of low-lying eigenvalues of the QCD Dirac operator in a finite 
volume, whose distribution can be predicted by Random Matrix Theory (RMT) 0, d E|- Many groups 0, S 0, 0, 
H, H, EH EH, El Ei, Ei, Ei, Q have computed S from fits to the distributions of low-lying eigenmodes of the Dirac 



operator. This seems to require less computer power for equivalent accuracy[17| than "asymptotic" methods, and also 
seems to involve fits to data with fewer free parameters. 

A method to compute F using eigenmodes was recently devised by G. Akemann, P. H. Damgaard, J. C. Osborn 
and K. Splittorff [18(. They imagine doing simulations in the epsilon-regime with two dynamical flavors, and then 
considering the spectrum of a partially quenched theory in which two non-dynamical flavors are coupled to an imag- 
inary isospin chemical potential. The correlator of the ordinary eigenvalues with those of the quenched ones involves 
the product of the isospin chemical potential and F. The required eigenmode distributions are easily computed from 
an ordinary dynamical ensemble. This paper is a pilot study of the technique at lattice spacings which are small 
enough that one can imagine quoting a continuum number. 

We can almost complete the calculation. The missing ingredient is actually analytic: We are doing the calculation 
in finite simulation volume. We know that finite volume alters the effective value of the condensate (said differently, 



the condensate appears in all relevant formulas multiplied by a volume-dependent coefficient) [l9j. We expect that F 
is subject to similar modification. The precise calculation of this connection does not appear to have been performed, 
and appears to require more knowledge of chiral perturbation theory than we have. We return to this point below. 



There have been two generations of previous variations of this idea, Refs. [20( and [2l|, |22j. These authors did 
simulations in which the isospin chemical potential was included in the action. This is of course much more expensive 
than what we are doing here. 
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We performed the calculation with lattice fermions which support the full SU(Nj) ® SU(Nf) chiral symmetry 
(including the anomalous singlet current and the index theorem), overlap fermions 23, 24j • We find doing overlap 
simulations to be reasonably straightforward. RMT makes predictions for QCD in fixed topological sectors and its 
predictions are for the pure imaginary eigenvalues of an anti-Hermitian Dirac operator. Both of these requirements are 
met unambiguously (modulo cutoff effects) by overlap fermions: the index theorem as realized by the Ginsparg- Wilson 
relation gives the topology and the usual definition of the condensate as £ = (1/ Z)dZ j 'dm stereographically maps 
overlap eigenmodes onto the imaginary axis 

Typical high-statistics supercomputer collaborations (for example Ref. quote f„ with 0.3 per cent errors (but 
slightly greater errors for F; see below). Ref. 2l|, which extracted F from eigenmodes (the isospin chemical potential 
is included in the simulation) has an 0.7 per cent error with about 10 3 samples. This was done deep in strong coupling 
so they did not attempt to make a continuum prediction. We will not come close to this amount of accuracy due to 
lack of statistics, but our results may not be uninteresting. Finding F does not seem necessarily to be a supercomputer 
project. 

If we were doing simulations in the real world of three light flavors, what would we expect? Two representative 
continuum chiral perturbation theory fits to data give F — 88 MeV [3] or F = 86.2 ± 0.5 MeV [291 ] . Our simulations 
are for two dynamical flavors and so strictly speaking we should not be talking about experimental comparisons. 
However, we will be unable to resist doing this. 

The rest of the paper is organized as follows: We first review the observables we are going to compute and the 
corresponding Random Matrix predictions. Then we discuss the systematic and statistical uncertainties of this 
simulation, focusing on the corrections due to the finite volume. After that, the setup and results of our actual 
simulation will be given. 



II. THEORETICAL BACKGROUND 

The key to the calculation is the connected correlation function between eigenvalues of the Dirac operator computed 
in the presence of a quenched chemical potential Aj, and ordinary eigenvalues computed at zero chemical potential, 
A,- 



(2) 



Eq. (3.49) of Ref. [l8j gives the RMT prediction for this correlator in terms of the usual rescaled variables x = 
Xi'SV, y = XjT,V, the rescaled mass m = m q Y,V, and the rescaled isospin chemical potential p, = 5 = fj,FyV. V is 
the volume. 



(2) conn / ~ ~\ ~ ~ 

P(i'i) {x,y)=xy 



l + (y,im) J v (im) imJ„+i(im) 
2 + (y,irh)' J v (irh)' (imJ„+i(im))' 
l + (y,x) J„(x) xJ v+ i{x) 



-I°{x,irh) J v (im) imJ v+ i(im) 
-l°{x,im)' J v [im)' {imJ v +\{im))' 



J v {im) imJ^ + i(im) 
J v -i{irri) imj v (im) 



(3) 



The prime indicates the derivative with respect to im. This formula is the limit of a more general one with two non 
degenerate quark masses, x and y are the rescaled eigenmodes for fi — and fi — S respectively. The integrals are 



I°(x,y) 



l ± {x,y) 



1 



dt J„(xVt)Jv(yVi) 



xJ v +i{x)J v {y) - yJ v+ x{y)J v [x) 



x 2 — y 2 



dte ±pt/2 J„(xVi)Jv(yVi) 



(4) 



Lattice data is never as beautiful as an analytic formula. In our case, the problems we face in applying Eq. [3] to 
simulation data are that 



• We have limited statistics. 
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FIG. 1: C(x) and I(X) from Eqs. Eland© for fj, = mEV = 2, ( max T,V = 7,v = 0. 



• We do not compute a large number of eigenmodes. 

• From our analysis of the integrated eigenvalue spectra we suspect that cutoff effects modify the spectrum of the 
higher modes away from the RMT prediction. 

• Simulations arc done in finite volume. 

Each of the four points contributes to the systematic or statistical error of our final result. 

The first problem on our list is the limited statistics. It makes the comparison of the measured eigenvalue data to 
the theoretical distribution of Eq. [3] difficult, ft is therefore convenient to compute integrated correlators: The one 
suggested by Refs. and [iH is 

C(x, Cmax) = / dyp(x +y,y) . (5) 



To avoid binning the data, which is highly ambiguous given the low statistics, we perform a second integral and take 
instead 

(max) = / C(X, (max)dx . (6) 

J-Cma* 

The first integral C(x) shows a spike at x — whose width goes roughly as 5 2 . Then I(X) will show a sharp step at 
X = 0. This is illustrated in Fig. [T] (In practice, we have only studied I(X).) 

To deal with the second and third problems, we realize that for any finite (max, only a finite number of eigenmodes 
contribute to C(x) or I(X). Since we will be measuring the eigenvalue correlations at tiny 8, a plot of the individual 
eigenvalue distributions will easily reveal which modes saturate the correlator below a certain value £max of A. An 
example is illustrated in Fig. [5] In this case mT,V — 2. We see that for (max'EV < 7 only the two lowest modes 
contribute to the v = correlator, and similarly for Qmax^V = 8 for \v\ = 1. 

Let us now come to the fourth problem. Finite volume typically modifies field theory correlation functions, because 
the propagators are altered by the boundary conditions. In the epsilon regime, finite volume also suppresses long- 
range order. To compute the correction to the condensate, one basically writes the Goldstone boson fields as U (x) = 
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FIG. 2: An illustration of the RMT eigenvalue density (for all modes) and the distributions of single eigenvalues, for /i = 
rriSV = 2 as a function of C, = AEV, in topological sectors (a) v — 0, (b) \u\ = 1. 



Co exp(in(x)/F), where Uo is the zero mode field which enters into all infinite- volume correlators (and is the variable 
connected to RMT formulas) and ir(x) is small (~ 1/L d+2 in d— dimensions). Then one expands the partition function 
to desired order in 1/L and integrates out the small terms. For the condensate, the effect is to replace E by Ej, = p^E 
where 

p s = l + C s -^A(0) + ... (7) 

with A(0) the contribution to the tadpole graph (propagator at zero separation) from finite-volume image terms. In 
the epsilon regime, A(0) = —D/\/V and D depends on the geometry [jof. (It is 0.1405 for hypercubes.) For Nf 
flavors, C* s = -(Nj -1)/N f . 

We are aware of no equivalent calculation appropriate to the case in hand. It amounts to computing corrections to 
derivatives of the partition function in the presence of a static background isospin vector potential. It includes both 
physical quarks (or pseudoscalar mesons) with ghost states. 

However, at the end of the day, we want to make a comparison with phenomenology. We will make a guess (readers 
who do not like plausibility arguments are invited to look away for two paragraphs) : 

The epsilon-regime correction to E is the same as the one- loop correction in the p-regime, with the same Cs. All 
that differentiates the two is the evaluation of A(0). We will just assume that the same result obtains for F - that 

F L = p F F (8) 

where pp is the analog of pp with Cf = — JV//2 replacing Cg. 

Is this a plausible guess? Epsilon regime axial vector and vector current correlators in boxes of length T have been 



computed by several authors 3l|. The mass- independent part of the correlator is a constant, 



F 2 

C(t) 



N f f D T 



1 + — o \-F= ~ k 

F 2 \VV V 



(9) 



In a symmetric box (V — T 4 ), k — D/2 and the effect is equivalent to Fl — PfF with Cf = —Nf/(2V2). Our guess 
scales properly with volume, but has a different factor. There is an additional 1/yV correction in C(i), which carries 
a nontrivial t dependence. It might be, that the effect of finite volume is to modify the functional dependence of Eq. 
|3l Clearly, we have exceeded the bounds of the numerical simulation we set out to do. 



5 



M 


N 


mode 1 


mode 2 


mode 3 





30 


0.34 


0.30 


0.31 


1 


75 


0.29 


0.33 


0.53 



TABLE I: The quality factor of the individual modes for the result of a combined RMT fit to the lowest mode in each topological 
sector, see Fig. [3] 



III. NUMERICAL SIMULATION 



Let us now come to the discussion of the actual simulation. Finding eigenmodes at nonzero isospin is easy: take a 
set of gauge configurations from the dynamical stream, rotate the time-like links 

U (x) -> exp(i[i)U (x) (10) 

and re-compute eigenmodes of the massless Dirac operator. This (standard [32| ) implementation insures that the 
chemical potential - and hence our extracted F - is not renormalized by lattice interactions. 

Our version of the overlap uses various implementations of fat links, but as Eq. [TU] is a global rotation, the fat 
links are rotated by the same amount. (That is, we first fatten the links, then rotate them.) And these quenched 
measurements of eigenvalues are basically "for free," as compared to the cost of the generation of the configurations. 

Our data set uses a lattice volume of 12 4 points. The lattice spacing a, as determined from the Sommer parameter 
tq [33} ] - is ro/a = 3.71(5). This calculation is "parasitic"; our main research program with our simulations involves 
heavier quark masses. We present data from the lightest quark mass at which we have extensive statistics, am q — 0.03. 
The pseudoscalar mass is am^ = 0.329(3). This is not in the epsilon regime, but experience shows RMT is robust 
enough to work well outside the epsilon regime, and anyway, we are just making a pilot study. In our conversions to 
physical units we use ro = 0.5 fm. 

The overlap operator uses a "kernel action" (the nonchiral action inserted in the usual overlap formula) with nearest 



and next-nearest (diagonal) neighbors. The 12 data set uses the differentiable hypercubic smeared link of Ref. [34 1. 



Details of the actions are described in Refs. 13|, [35j, [3fj, [37[; everything except the choice of smearing for the link is 
identical to previous work. 

Let us briefly summarize the algorithmic set-up of our computation. For the Nf — 2 simulations we use the 
reflection/refraction algorithm first devised in Ref. [3c| . In order to improve the tunneling rate and precondition 



the fermion determinant we use two additional heavy pseudo-fermion fields as suggested by Hasenbusch 39( . The 
integration is done with multiple-time scales j4jj. The runs for all sea quark masses were performed within a few 
months on a cluster of 32 Opteron CPU's which are connected by an Infiniband network. For the am q = 0.03 
ensemble, we collected about 400 thermalized HMC trajectories of unit length and analyzed 30 v — lattices and 75 



\v\ = 1 ones. We are computing eigenmodes using the "Primme" package of McCombs and Stathopolous 41 1. 

Our analysis begins with a determination of a 3 £z, by fitting the RMT prediction of the distribution of the lowest 
eigenvalues of the Dirac operator without chemical potential to the data. As mentioned in the introduction, this 
method has been used by us and other groups during the last years and will serve as a cross check for the new 
method. Specifically, we fit the integrated eigenvalue density (cumulant) 

P<(A) = / Pi(X')d\' (11) 







with Pi(X) the density of the i-th lowest eigenmode. An example of a typical fit, this one to the lowest eigenmodc 
in the v = and \v\ = 1 sectors, is shown in Fig. [3] Its best fit value is cl 3 Y,l = 0.0073(3), the quality factors for 
each individual mode from the combined fit are listed in Table HI Fits to other sets of eigenmodes produce entirely 
consistent results. 

Now we proceed to the analysis of the eigenvalues with chemical potential and fits to Eq. [5J The correlator depends 
on two free parameters, 5 and the scaled condensate (times volume) T,V. The "experimental" I(X) varies almost 
continuously in X when there are many lattices in the data set. Since it is an integral, its values at different X 
are highly correlated. We account for correlations via a bootstrap analysis: We pick a bin size for X and make a 
set of bootstrap-averaged determinations of I(Xj, Cmax), with an uncertainty of a bin value Oj determined from the 
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0.00 0.02 0.04 0.06 0.08 0.10 



FIG. 3: Cumulant distribution for the three lowest eigenmodes (without chemical potential) in the \nu\ = and 1 sectors, for 
our 12 4 [3 = 7.3 am q — 0.03 data set. The curve is a fit to the lowest mode in v — and in \v\ — 1, with a 3 Ez, = 0.0073. 



bootstrap. The we do a fit minimizing 

We determine uncertainties in the fit parameters from the fluctuations in a bootstrap average over a number of these 
fits. The fits involve a number of arbitrary choices, and our results should be independent of them. We have tested 
all of the following 

• Range of ( max 

• Range of X. 

• Number of values of Xj 

The fit then determines a favored value of 6, from which Fl = 5/(VVfi). It also gives a second determination of 
the condensate T,lV, to be compared to the previously discussed direct fit to cumulant distributions. 

The cumulant distribution of the eigenmodes tells us how many modes contribute for any ( max - In this case, we 
believe that we can keep up to three eigenmodes and be consistent with random matrix theory predictions, so a 
maximum value will be someplace near the average value of the third eigenmode. This corresponds to a (not rescaled) 
(max = 0.07 for v = and 0.08 for \v\ = 1. 
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FIG. 4: Data (squares) and fit from \v\ = 1. A cut Cmax = 0.07 is enforced. The best fit values for this bootstrap sample are 
E L V = 147 and aF L = 0.071. 
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FIG. 5: Histogram of 200 bootstrap fits to aFL from the \v\ — 1 data set: f„ 



: 0.08 and X = 0.004 in the fits 



Fig. [6] shows the dependence of fits on Cmax- These fits all include only the three lowest (but nonzero) eigenvalues in 
a topological sector. Smaller Cmax means that less data is included in the fit; larger Cmax means that more eigenmodes 
are needed to saturate the correlation function. It appears that results for this data set do not depend on the choice 

of Cmax • 

Finally, a histogram of bootstrap values of the finite- volume decay constant Ft is shown in Fig. [5l 

As the best-fit values do not show strong dependence on Cmax we choose to quote results from the Cmax — 0.08 fit 
for the \v\ = 1 data set and Cmax = 0.07 for v = 0. They give aF L = 0.065(5) and a 3 S L = 0.0068(2); the \v\ = 1 set 
gives aF L = 0.073(5) and a 3 E L = 0.0072(2). Combining these fits gives aF L = 0.069(4) and a 3 E L = 0.0071(1). The 
fit gives a corroborating result for lattice-regulated S, consistent with our result from eigenmode cumulants. 

To complete the calculation of S, we need to convert from lattice regularization to MS. We do this using the 
intermediary of the RI-MOM scheme [3]. Our methodology is basically identical to what we did in Ref. [i~3l |. 

We collected data from the three quark masses, am q = 0.10, 0.05, and 0.03. We took 40, 30, and 40 lattices, 
respectively. In RI-MOM, a Z-factor is basically the product of an averaged unamputated momentum-space vertex 
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FIG. 6: Best-fit values of oFl, varying („ 



function times two averaged inverse propagators, so for an N-lattice data set we formed N single-elimination jackknifed 
vertex-propagator combinations. The jackknife-averaged Z-factors are shown in Fig. [7Ja). 

We want the Z-factor at a fiducial scale /i which conventionally is chosen to be 2 GeV. We fit the Z-factor for a 
particular mass to a linear dependence in [i for /it's "nearby" /i = 2 GeV. There is one value for each possible value 
of lattice momentum. Combining all the points with the same value of J^- kf gives a number with an error bar for a 
given value of J2i f° r each of these jackknifed data sets. 

We interpolate these fits to the desired value of /z, 2 GeV with the lattice spacing taken from our Sommer parameter 
measurements. The uncertainty comes from jackknifing these results. For our simulations, — 2 GeV corresponds to 
afi ~ 1.4, and we did fits over the range 1.0 or 1.1 to 1.6. (To be precise, 1.497 at am q = 0.10, 1.45 at 0.05, 1.368 at 
0.03.) The interpolated results are shown in Fig. [7Kb). A linear fit gives the am q = Z-factor of 0.656(28). 

To convert to MS scheme, we use the coupling constant from the so-called "ay" scheme. The one-loop expression 
relating the plaquette to the coupling is 



bx^TrU p = -^-a v (q*)W, 



(13) 



where W = 0.366 and q*a = 3.32 for the tree-level Liischer-Weisz action. All three data sets gave ay{q*) = 0.188. 
The lattice spacing does not vary much in the three data sets, so converting any of them gives a^ 1 (2 GeV)= 0.207 
and Z§ rs /Z§ 1 ' = 1.158. Thus the combined lattice to WS conversion factor is Zf s {2 GeV)=0.76(3). We insert our 
lattice determination for aF^ in p^. Dividing it out, we find 



r Y,(MS,(j, = 2 GeV) 1/3 = 0.594(13) 



(14) 



(S(M5,|i= 2 GeV)) 1/3 = 234(4)MeV 



(15) 



From the cumulant analysis, r £(MS,,u = 2 GeV) 1 / 3 = 0.600(13)or JZ(MS, fi = 2 GeV)) 1 / 3 = 237(6) MeV. This 
seems to be quite consistent with other recent determinations 3, 14, 15, 17 1. 
Now for F. Combining results for v = and \v\ — 1, we find 



r a F L = 0.255(13) 



(16) 



or with rn = 0.5 fm, 



F L = 101(6)MeV 



(17) 
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FIG. 7: (a) Zs from our RI-MOM analysis: octagons are am q — 0.10, squares, 0.05 and diamonds 0.03. Errors are from a 
jackknife average, (b) Zs(RI') interpolated to pi = 2 GeV, vs am q . The diamond shows the extrapolated point (with its error). 



This is a little high compared to the phenomenological estimates of Refs. [28|, 129J. There are ample reasons to 
suspect that some of this discrepancy could be due to fundamental problems with the simulation (Nf = 2 is not real 
world), but we have already argued that we expect that Fl differs from F, by a factor which differs from unity by 
0(1/(F 2 L 2 ). 

Using the guess for finite value corrections, Fl — plF, and our guess for pl, we estimate 

r F = 0.213(11) (18) 

or 

F = 84(5)MeV. (19) 

Perhaps the error is more interesting than the central value: We have produced a six per cent statistical uncertainty 
for a full QCD quantity, with an extremely modest investment of computer power. 

We remarked earlier that large scale simulations using asymptotic correlators typically produce calculations of / w 
with small uncertainties but noisier predictions for F. As an example, MILc[3| quotes (their Table IV, Fit A, our 
conversion) r\F = 0.131(10) or F = 82(6) MeV. {r\ is defined through the heavy-quark force as rfF(r\) = 1.0.) 
ETMC gets tqF — 0.197(2) [iH]. We believe that our result used several orders of magnitude less computer power 
to produce a comparably-sized statistical error. Control over the systematic error due to the finite volume, however, 
remains a topic of future research. 



IV. CONCLUSION 



This was a pilot study, and it needs several ingredients before a complete calculation of F can be performed: 

• A real analytic calculation of the finite volume correction to the eigenvalue correlator is needed. Is it of the form 
Fl = plF as we assumed, and if so, what is Pf? In the absence of such a calculation, simulations at several 
volumes will be needed. 

• A real numerical calculation needs to truly enter the epsilon regime, through a combination of lower pseudoscalar 
mass and larger volume. 

• The best way we know to extract the condensate from Dirac eigenmodes used the distribution of individual 
eigenmodes, not the sum over many modes. The analog for F would be a calculation of the correlator of 
one eigenmode of ordinary fermions with one eigenmode of partially quenched fermions in an isospin chemical 
potential. We understand [46j that such a calculation is in progress. 
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• Finally, while we said at the beginning of this paper that F apparently does not need to be a supercomputing 
project, finding the physical decay constants involves determining higher-order coefficients in the chiral La- 
grangian. At present, this can only be done in a spectroscopy-type calculation, from the asymptotic behavior 
of a correlation function. This can be expensive. Is there a way 47| to remove that bottleneck by relating these 
coefficients to some property of Dirac eigenmodes? 
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